#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "MSTK.h"
#include <pthread.h>
#include <assert.h>

#ifdef __cplusplus
extern "C" {
#endif

double tictoc (int n);


  int MESH_Smoothing_Laplacian(Mesh_ptr mesh, int weight, int boundary) {
  int i, j, idx;
  int  nvv;
  int nv, ne, nf, nr;
  MVertex_ptr mv, vertex;
  MEdge_ptr me;
  MFace_ptr mf, mrfaces;
  MRegion_ptr mr;
  List_ptr mfverts,mrverts,mvertices;;
  
  nv = MESH_Num_Vertices(mesh);
  ne = MESH_Num_Edges(mesh);
  nf = MESH_Num_Faces(mesh);
  nr = MESH_Num_Regions(mesh);
 
  int element_dim = -1, boundary_dim = -1;
  int num_element = 0;
  if (nr) {
    element_dim = 3;  boundary_dim = 2;
    num_element = nr;
  }
  else if(nf) {
    element_dim = 2; boundary_dim = 1;
    num_element = nf;
    }
  else {
    fprintf(stdout,"\nThis is not a valid mstk file for Laplacian smoothing\n");
    exit(-1);
  }  
      /* decide if it is 2d planary mesh or 3D surface mesh */
    if (element_dim == 2) {
      double xyz_coord[3];
      double xyz_coord0[3];
      vertex = MESH_Vertex(mesh,0);
      MV_Coords(vertex,xyz_coord0);
      for (j = 1; j < nv; j++) {
	vertex = MESH_Vertex(mesh,j);
	MV_Coords(vertex,xyz_coord);
	/* if z coordinate is not all the same, it is a surface mesh */
	if ((xyz_coord[2]!=xyz_coord0[2]))
	  boundary_dim = 2;
      }
      if (boundary_dim == 1)
	fprintf(stdout,"\nThis is a 2D planary mesh\n");
      
    }

  /* Loop through vertex */
  idx = 0;
  double coords[3],vcoords[3];
  tictoc(1);
  while((mv = MESH_Next_Vertex(mesh,&idx))) {
    if(MEnt_GEntDim(mv)>boundary_dim) {
      mvertices = MV_AdjVertices(mv);
      nvv = List_Num_Entries(mvertices);
      MV_Coords(mv, vcoords);
      /*
      fprintf(stdout,"vertex %d has %d neighbors\n",MEnt_ID(mv),nvv);
      fprintf(stdout,"before: vertex %d (%lf,%lf,%lf)\n",MEnt_ID(mv),vcoords[0],vcoords[1],vcoords[2]);
      */
      coords[0]=0.0;
      coords[1]=0.0;
      coords[2]=0.0;
      for (i = 0; i < nvv; i++) {
	vertex = List_Entry(mvertices,i);
	MV_Coords(vertex, vcoords);
	coords[0] += vcoords[0];
	coords[1] += vcoords[1];
	coords[2] += vcoords[2];
      }
      coords[0]/=nvv;
      coords[1]/=nvv;
      coords[2]/=nvv;
      MV_Set_Coords(mv,coords);
      MV_Coords(mv,vcoords);
      /* 
     fprintf(stdout,"after : vertex %d (%lf,%lf,%lf)\n",MEnt_ID(mv),vcoords[0],vcoords[1],vcoords[2]);
      */
    }
  }
  fprintf(stdout,"time is %lf \n", tictoc(2));
  return 1;
  }

  

#ifdef __cplusplus
}
#endif

